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ABSTRACT 

We have developed a new galactic chemo-dynamical evolution code, called GCD+, 
for studies of galaxy formation and evolution. This code is based on our original three- 
dimensional tree N-body/smoothed particle hydrodynamics code which includes self- 
gravity, hydrodynamics, radiative cooling, star formation, supernova feedback, and 
metal enrichment. GCD+ includes a new Type II (SNe II) and la (SNe la) supernovae 
model taking into account the lifetime of progenitor stars, and chemical enrichment 
from intermediate mass stars. We apply GCD+ to simulations of elliptical galaxy for- 
mation, and examine the colour-magnitude relation (CMR), the Kormendy relation, 
and the [Mg/Fe]-magnitude relation of simulation end-products. GCD+ is a useful and 
unique tool which enables us to compare simulation results with the observational data 
directly and quantitatively. Our simulation confirm the results of Kawata (2001) who 
uses a simpler chemo-dynamical evolution code. We newly find that radiative cooling 
becomes more efficient and thus the gas infall rate increases, with decreasing mass of 
galaxies, which contributes to the slope of the CMR. In addition, the sophisticated 
treatments of both SNe II and SNe la in GCD+ show that feedback from SNe la plays a 
crucial role in the evolution of elliptical galaxies. We conclude that the feedback effect 
of SNe la should not be ignored in studying the evolution of elliptical galaxies. 

Key words: galaxies: elliptical and lenticular, cD — galaxies: formation — galaxies: 
evolution — galaxies: stellar content 



INTRODUCTION 
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to place within galaxies. Since heavy elements are believed 
to have been synthesized in stars, the metal abundances in 
present-day galaxies should offer a record through which we 
may trace the history of star forma t ion w ithin galaxies. Since 
the pioneering work of Tinsley (1972), studies of chemi- 



cal evolution have succeeded in connecting observed metal 
abu ndances to the history of star formation within galaxies 
e.g. lArimoto fc YoshTi] |l987|; [Matteucci fc Tornambc ||l987t 



Timmes, Wooslcy, fc Weaver | |1995|; [Gibson 1 11997 



pini, Matteucci, fc Gratton 1997; Tantalo et al 



fe1|T98 

HChai 
L998b 



Moreover, recently developed semi -analytic models (White 



h Rees 



1987 



White fc Frenk 1991) including chemical evo- 



lution make it possible to relate the observed chemical prop- 
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20011: Somerville, Primack, fc Fabcr 2001; [Beasley et al 
^1 



20021). These techniques require low computational costs 



and allow large coverage of parameter space. Although the 
pure chemical evolution models and the semi-analytic mod- 
els seem to succeed in explaining the observational proper- 
ties of galaxies at various redshifts, they inevitably assume 
a phenomenological model involving a number of parame- 
ters to describe the processes of gas accretion rates, radia- 
tive cooling, star formation, and supernova feedback within 
a galactic halo. In addition, structures of the end-products 
can not be discussed, because there is no information about 
the dynamical history. 

On the other hand, recent advances in computer tech- 
nology and numerical methods have made it possible to cal- 



erties of galaxies to cosmological theory (e.g. 


Kauffmann fc 


consis 








2001; 



2001; [Chiosi fc Carraro" 2002). Three dimensional chemo- 



dynamical evolution codes have been applied to the foUow- 
ing studies over the past five years: Steinmetz fc Mxiller 
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(1995) studied the chemo-dynamical evolution of disk galax- 
ies and succeeded in distinguishing the che mical properties 



Chiba 



Berczik 



between hal o, bulge, and disk stars (see a l sopekki 
200C|, ^00 1|) . [Raiteri, Villata, fc Navarro | dlOgel ) and 
(|1999|) took into account metal enrichment by Type la su- 
pernovae (SNe la) as well as Type II supernovae (SNe II), 
and reproduced the correlation between [O/Fe] and [Fe/H] 
for stars in the solar neighbourhood. Such self-consistent 
calculation of chemical evolution allows one to analyse the 
photometric properties of simulation end-products in combi- 
nation with stellar population synthesis, and to obtain the 
luminosity of the end-products without the assumption of 



metz 



ique, steinmetz & Navarro (1999), 


Navarro & Stein- 


(200C), and Koda, Sofue, & Wada 


(I2000D discussed 



the zero-point of the TuUy - Fishe r relation o f disk galax- 



Bekki fc Shioya (1998, 



1999 



2000 



2001) studied the 



dynamical and photometric properties of elli ptical galax- 
ies formed by a merger of two disk galaxies. Mori et al. 



(1997) showed that the galactic wind caused by supernovae 
(SNe) feedback leads to a positive colour gradient observed 
in so me dwarf ellipt ic al galaxies. Using cosm ol ogical simula- 



tioHj^ IGnedin 



(20011), and Theuns et al. 



(tl998|), ICen fc Ostriker 



(1999), Tissera et al. 



of the intcrgalactic medium (see also Aguirre et al 



(2002) st udied metal enrichm ent 

as 



2001 



a w ork coupling semi-analyt ic model to hydro simulations). 



and 



Nagamine et al. 



(2001) examined the luminosity func- 
tion and colour distribution function of galaxies at different 
redshifts. Since radiative cooling depends on the metallic- 
ity, dyna mical evolution is quite sensiti v e to the chemica l 

poool). 



evolution (Kaellander fc Hultman 



1998 



Kay et al. 



Thus, we must calculate chemical and dynamical evolution 

self-consistently. 

Us ing a chemo-dynamical evolution code, Kawata 
( ^001b| , hereafter KOlb) were able to reproduce the colour- 



magnitude relation (CMR) of observed elliptical galaxies. 
KOlb showed that SNe feedback affects the evolution of 
lower mass systems more strongly, and induces a stronger 
galactic wind in lower mass systems. Such low mass sys- 
tems have lower metallicities and bluer colours than higher 
mass s ystems, as pr e dicted by ana lytic models (e.g. 



moto & Yoshii 



1987; 



Gibson 



Ari- 



1997 ) and more phenomeno- 
[Carlberg 



1974 



logic al numerical simulations (e.g. Larson 
1984). Moreover, KOlb found that galactic winds are trig- 
gered mainly by SNe la rather than SNe II, although KOlb 
ignored the lifetime of progenitors of SNe II and SNe la 
(i.e. KOlb assumes the instantaneous recycling approxima- 
tion for SNe II, and all SNe la occur simultaneously 1.5 Gyr 
after stars were born). 

We have developed a new chemo-dynamical evolution 
code, called GCD+, which adopts a more sophisticated chem- 
ical evolution model than that discussed in KOlb. We have 
relaxed the instantaneous recycling approximation for SNe 



II, and adopt a more sophisticated SNe la model (Kobayashi 



Tsujim jto, & Nomoto 2000, hereafter KTNOO). Our code 



t akes into accou nt the metallicity dependent lifetime of stars 
( |Kodama||l99^), and metal enrichment from intermediate 
mass stars ( [van den Hoek fc Groenewegen 1997). The pur- 



pose of this paper is to introduce details of GCD+. In addition, 
as an application of GCD+, we carry out simulations of ellip- 
tical galaxy formation similar to the ones of KOlb. Then 
we re-examine the CMR, the Kormendy relation, and the 
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Figure 1. Cooling rates as a function of temperature. Each 
curve corresponds to a cooling rate with different metallicities as 
indicated. 



[Mg/Fe]-magnitude relation, which were studied in KOlb. 
Since these observed relations are well-established, and are 
expected to reflect the chemo-dynamical evolution of ellipti- 
cal galaxies, this is an important application for GCD+. Next, 
taking advantage of the sophisticated SNe II and SNe la 
models in GCD+, we clarify the role of SNe la, compared to 
SNe II, in the evolution of elliptical galaxies. 

The outline of this paper is as follows. Details of GCD+ 
are described in Section ^ In Section ^ we use GCD+ to sim- 
ulate elliptical galaxy formation. We show the galaxy for- 
mation model used in this paper briefly in Section 3.1. In 



Sections 3.2 



-3.6, results of numerical simulations are pre- 
sented, and are compared to the observed scaling relations. 
Discussion about the study of elliptical galaxy formation is 
given in Section Finally, we present our conclusions from 
this paper in Section ^ 



THE NEW CHEMO-DYNAMICAL 
EVOLUTION CODE 



Kawata 


(1999) and KOlb. The code is essentially based on 


TreeSPH (Hcrnquist fc Katz 198^; Katz, Weinberg, fc Hern- 


quist 


19961), which combines the tree algorithm (Barnes & 


Hutt 


198h) for the computation of the gravitational forces 



with the smoothed particle hydrod ynamics (SPH: Lucy 



1977; Gingold fc Mnoraghan 1977) approach to numerl 



cal hydrodynamics. The dynamics of the dark matter and 
stars is calculated by the N-body scheme, and the gas com- 
ponent is modeled using the SPH. It is fully Lagrangian, 
three-dimensional, and highly adaptive in space and time 
owing to individual smoothing lengths and individual time 
steps. Moreover, it includes self-consistently almost all the 
important physical processes in galaxy formation, such as 
self-gravity, hydrodynamics, radiative cooling, star forma- 
tion, SNe feedback, and metal enrichment. The code is vec- 
torized and parallelized. The parallelization is done using 
the MPI library, so that the code can run on any type of 
computer, including supercomputers and PC clusters. We 
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have revised the modehng of cooUng, star formation, and 
SNe feedback drastically from the code used in KOlb. Here, 
we will only deal with the new aspects of our N-bo d y/SP H 
code, which has already been described in 



LO 



Kawata 



( |1999| ). 



2.1 Gas Cooling 

To model gas cooling, we use the cooling curves computed by 
MAPPINGS IIl[] written by R.S. Sutherland. Using MAP- 
PINGS III, we have computed the cooling function of the 
ionization equilibrium gas with various metallicities. Fig. ^ 
shows all the cooling curves which we use. With linear in- 
terpolation of these cooling curves, the code calculates the 
cooling rate of A(T, [Fe/H]) as a function of the tempera- 
ture and metallicity for each gas particle. We assume that 
the cooling rate of the gas with metallicity of lower than 
[Fe/H] = —3 is the same as the one with [Fe/H] — —3, and 
the cooling rate of the gas with [Fe/H] > 0.5 is equivalent 
to that at [Fe/H] = 0.5 to avoid extrapolation. The cooling 
rate for the gas with the solar metallicity is larger than that 
for the gas of primordial composition by more than an order 
of magnitude, thus cooling by metals should not be ignored 
in numerical simulations of galax y formation (Kaellander & 



Hultman 



199^ 



Kay et al. 2000). The temperature for each 
Pi^j,mp/{kBpi), where Pi and pi, 
i-th particle, and fi, ku, 
Boltzmann's con- 



particle is derived by Ti 
are the pressure and density for 
and nip are the mean molecular weight 
stant, and the proton mass, respectively. For simplicity, we 
fix p = 0.6, regardless of the metallicity. We set the lower 
limit of the temperature to be Tlim ~ 100 K. 



2.2 Star Formation and Initial Mass Function 

We model st ar for r natio n usi ng a method similar to that 



Katz 



( |199: 



and 



Katz, Weinberg, fc Hernquist 



s ugges ted by 

(1996). The criteria for star formation and the star formation 
rate are the same as those in KOlb. We use the following 
three criteria for star formation: 1) the gas density is greater 
than a critical density, pcrit = 2 x 10"'^^ g cm"'', i.e. nu ~ 
O.lcm"^; 2) the gas velocity field is convergent, V • < 0; 
and 3) the Jeans unstable condition, h/cs > tg, is satisfied, 
here h, Cs, and tg = y^Sir/lGGpg are the SPH smoothing 
length, the s ound speed, ari d the dynamical time of the gas 
respectively (Kawata 1999). When a gas particle is eligible 
to form stars, its star-formation rate (SFR) is 



dp, 
~dt 



dps 
dt 



c*Pe 



(1) 



where c* is a dimensionless SFR parameter and tg is the 
dynamical time, which is longer than the cooling timescale 
in the region eligible to form stars. This formula corresponds 
to the Schmidt law that SFR is proportional to Pg'^. We fix 
c, = 0.5, following KOlb. 

We assume that stars, which are repr esented by a sta r 
particle, are distributed according to the salpeter (195 
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scrihcd in |Siithcrlarid Donita. _il22^1^and available at 



http: / /www. mso.anu.edu.au/~ralph/map. html 
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Figure 2. Chemical yields as a function of age of a star particle 
with mass of 1 Mq. The upper (lower) panel shows the total 
ejected oxygen (iron) mass. The thick (thin) line indicates the 
history of a star particle with the metallicity of Z = 0.02 (10~*). 



initial mass function (IMF). The IMF by number, "I'(m), in 
each mass interval dm is defined as^ 



${m)dm = Am 



dm, 



(2) 



where x = 1.35 is the Salpeter index and the coefficient 
A is determined by the normalization in the mass range 
Ml <m < Afu. We set Mi = 0.2 Mq and AU = 60 Mq in 
this paper. 



2.3 Feedback 

GCD+ takes account of energy feedback and metal enrich- 
ment by SNe, and metal enrichment from intermediate mass 
stars. We consider here both SNe II and SNe la. The code 
calculates the event rates of SNe II and SNe la, and the 
yields of SNe II, SNe la, and intermediate mass stars for 
each star particle at every time step, considering the IMF 
and stellar lifetimes. Here, we assume the same stellar life- 
times as the ones used in 



Arimoto 



Kodama 



(1997) and Kodama & 



(1997). The simulation follows the evolution of the 



Equation (H) corresponds to an IMF by mass oc m 



4 D. Kawata and B.K. Gibson 



abundances of several chemical elements C^H, ^Hej^^C, ^''N, 
^•^0, 2°Ne, 2*Mg, 2*Si, ^"^Fe, and Z, w here Z means the to 



t al me tallicity) . We refer to Table 2 of Woosley 
(1995) for the solar abundances. 



Weaver 



2. 3. 1 Type II Supernovae 

We assume that each massive star (> 8 Mq) explodes as 
a Type II supernova. To calculate the ejected mass of gas 
and metals by SNe II, we use the stellar yields derived by 
Woosley & Weaver (1995). Woosley & Weaver (1995) pro- 



vide yields for a grid of stellar masses between 11 and 40 
Mq, and metallicities between and 1 Zq. With linear in- 
terpolation of these grids, we ob tain vields as a func t ion o f 



the stellar mass and metallicity. Woosley fc Weaver (199. 



considered three different models for m > 30 Mq, their A, 
B, and C models. These differ in the amount of energy im- 
parted by the piston in their models at explosio n initiation. 
Following Timmes, Woosley, fc Weaver ( 1995 ), we use the 
'B' model in this mass regime. For m > 40 Mq, we as- 
sume the same abundance ratios for the yields and the same 
mass fraction for the remnants as those for 40 Mq star. For 
Z > Zq, we use the following simple scaling, e.g. the yield 
of carbon from a star with the mass, m, and metallicity, Z, 



mej,c{Z,m) = 



(3) 



mej,c{ZQ,m) +mej{ZQ,m)Zc^(i 
x{Z/Zq-1), 

where mej,c is the total mass of ejected carbon, 
ing newly synthesized and initial carbon; rriej — ? 
rriej^He + mej.z indicates the total ejected mass; Z, 
the solar abundance of carbon. The yields from stars with 
initial masses between 8 and 11 Mq ar e still unclear (e.g. 

Hence, we simply 



includ- 

1ej,H + 



c,0 



is 



Hashimoto, Iwamoto, & Nomoto 



1993) 



interpolate linearly b etween the yields for the lo west mass, 
i.e. ~ 11 Mq, star in Woosley & Weaver (199E ) and those 



for the highes t mas s, i.e. 8 Mq star in van den Hoek & 



Groenewegen ( 1997| ), as mentioned in the next section 



2.3.2 Intermediate Mass Stars 



The yields and remnant m asses for intermediate mass stars 
( < 8 M q) are taken from [van den Hoek fc Groenewegen 



( tl997D 
yields of 



van den Hoek & Groenewegen 



(1997) present the 



^H, *He, ''^C, "C, "''N, and ^^OTor stars with 
initial masses between 0.8 and 8 Mq and initial metallic- 



ities of Z = 0.001 - 0.04. We d o not use "C yields, van 



den Hoek & Groenewegen ( 1997 ) provide the n e wly fo rmed 
and ejected metals unlike Woosley & Weaver (1995) who 



present the total ejected metals. The refore, we calculate the 
total eject e d me tals using Table 1 of van den Hoek & Groe- 
newegen (1997) for H, He, C, N, and O and the initial 
abundances for the other elements, i.e. Ne, Mg, Si, and Fe, 
which are simply scaled to the solar abundance set. As for 
stars with Z > 0.04, we use the same scaling as equation 
(|) based on ^ = 0.04 yields. For Z < 0.001, we simply as- 
sume the same yields as ^ = 0.001 yields. Whilst this sim- 
ple assumption may overestimate the ejected metals, since 
the metals are ejected from intermediate mass stars later 
than those from SNe II, which eject much larger amounts of 
metals, this simplification should not affect the final metal- 
licity. Nevertheless, the knowledge of the yields for low and 



zero metallicity stars is important in studies of the chem- 
ical composition of extremely metal poor objects observed 
at both low and hi gh redshifts, although i t is n o t yet well- 



al. 



2001 



establishe d (e.g. Fujimoto, Ikeda, fc Iben 



Marigo et al. 2001) 



200C 



Chieffi, et 



2.3.3 Type la Supernovae 

We adopt the SNe la model proposed by KTNOO, which 
is based on the SNe la progenit or scenario proposed by 



Hachisu, Kato, & Nomoto ( 199E| ). This model is difi'crcnt 



f rom t he conventional one proposed by Greggio & Renzini 



( [1983| ). Details of this model are described in KTNOO. Here, 



we briefly explain this model and implementation in our 
code. Following KTNOO, SNe la are assumed to occur at bi- 
nary systems which consist of primary and companion stars 
with appropriate masses and metallicities ([Fe/H] > —1.1). 
The primary stars have the main-sequence mass in the range 
of irip^i = 3 Mq and nip^u = 8 Mq , and evolve into the C-l-O 
white dwarf (WD). The mass range of companion stars are 
restricted to between nid.RGA = 0.9 Mq and md,RG,u = 1.5 



Mq and between md,_ 



= 1.8 Mq and md,_ 



= 2.6 

Mq. KTNOO designate the binary system whose compan- 
ion has the mass within the former (latter) mass range 
"RG-f-WD (MS-I-WD) system". Finally, the total number of 
SNe la is obtained by the following equation as a function 
of age, t, of a star particle with the mass of rris M, 



TVs 



dm 



max(m^_MS,!>™t) 



0' 



m ^dm 



Ml 

^d{'m)dm 



pd «G,„ ^Jm)dm 



(4) 



where mt is the mass of star whose lifetime is equal to t, and 
the mass function of the companion stars is assumed to be 
^d{m) DC m"^-^^ by number, following KTNOO. The term be- 
fore the square bracket indicates the number of C-l-O WDs, 
i.e. primary stars. The first term within the square bracket 
determines how much fraction of C-l-O WDs evolves into 
SNe la from the MS-f WD systems, and the second term de- 
termines how much fraction of C4-0 WDs evolves into SNe 
la from the RG-I-WD systems. Following KTNOO, we set 
feMS = 0.05 and &rg = 0.02. The nucleosynthesis prescrip- 
tions for SNe la are taken from the W7 model of Iwamoto 



et al. 



(|1999[). 

Fig. ^ shows the total amount of oxygen and iron ejected 
from a star particle with the mass of 1 Mq as a function 
of its age. Initially, metals are ejected only by SNe II. This 
continues until the 8 Mq star dies (~ 0.04 Gyr in the case 
oi Z = 0.02). Oxygen is produced mainly by SNe II. After 
SNe II cease, the continuous ejection of oxygen is mainly 
due to the contribution from intermediate mass stars. On 
the other hand, a significant fraction of iron is produced by 
SNe la. Stars with Z = 0.02 lead to SNe la depending on the 
lifetime of their companion stars. SNe la occur in MS-I-WD 
(RG-I-WD) systems in the range of age of 0.7 - 1.5 (> 2.8) 
Gyr. However, a star particle with Z = 10"'* does not lead 
to SNe la, because the above SNe la model restricts the 
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metallicity range for progenitors of SNe la to [Fe/H] > —1.1. 
For simplicity, we assume that the metaUicity range for SNe 
la is logZ/ZQ > -1.1, instead of [Fe/H] > -1.1. In the 
implementation, the code possesses a look-up table of the 
yields of all the chemical elements, remnant masses, number 
of SNe as a function of the age and metallicity, and calculates 
those values for each star particle at every time step. 



2.3.4 Energy Feedback 

One of the most difficult and most critical processes to model 
in galaxy formation simulations is the way in which the en- 
ergy feedback from SNe affects the surrounding gas. Unfor- 
tunately, there is no clear understanding of how thi s should 
be modeled . We a dopt a simple model proposed by Navarre 
(|l993[) 



fc White (1993) and used in KOlb. This model assumes 
that the energy produced by SNe affects only the tempera- 
ture and velocity field of the surrounding gas, and its effect 
is implemented by increasing the thermal (Eth) and kinetic 
(^'kin) energy of the gas neighbours of each star particle by 
an amount corresponding to the energy released by SNe. We 
assume that a parameter /„ = -E'kin/(£^kin + Sth) defines the 
fraction of the available energy to perturb the gas velocity 
field, and the rest of the energy of the SNe contributes to the 
increase in the thermal energy of the gas (see KOlb for de- 
tails). It is known that kinetic feedback affects the history of 
star formation more strongly than thermal feedback, which 
quickly dissipates due to efficient radiative cooling where the 
gas density is high enough to form stars. The parameter, /t,. 



controls the ma gnitude of the effect of SNe (KOlb; Navarre 



fc White 1993). We assume that each SN yields the energy 



of esN X 10 ergs, and then i?kin ~ fvf^SNW ergs. Because 
an initial SN energy has not been established quantitatively 
yet, we consider the available SN energy to be a free param- 
eter. Finally, in our code, there are two parameters, and 
esN, to control the magnitude of the effect of SNe. 



2.3.5 Implementation 

Based on the above feedback processes, the code calculates 
the amount of the mass, energy, and heavy elements re- 
leased from each star particle within each time step. Our 
c ode a dopts the simple feedback scheme suggeste d by [Katz 



; ode a dopts the simple teed back scheme suggeste d by Kat; 
(1992) (see also Lia, Portinari, & Carraro ^002| , for a pro- 



posed alternative scheme). The mass, energy, and heavy el- 
ements are smoothed over the neighbouring gas particles 
using the SPH smoothing algorithm. For example, when the 
z-th star particle ejects the mass of MsN,i, the increment of 
the mass of the j-th neighbour gas particle is given by 



AMsN,- 



where 



-MsN,,W{rij/h^), 



and W(x) is an SPH kernel. 



8 



W^.) = ^A 2(1-.) 



1-6x^-^62;^ if0<a;<l/2, 
^ ifl/2<a;<l, 
otherwise. 



(5) 



(6) 



(7) 



Here, Tij is the distance between the i-th and j-th particles, 
and hi is the smoothing length for the i-th particle. Note 
that the above equations employ hi, instead of hij = {hi + 
hj)/2 which is adopted in other S PH calculations , such as 
density and pressure gradient (see Kawata 1999). This is 



because feedback is a one-way process, which does not have 
to be symmetrised. In addition, if all neighbour particles 
have much smaller smoothing length than the i-th particle, 
hij is possible to become smaller than rij , and then the i-th 
particle loses the place to feedback. The use of hi ensures to 
avoid this problem. 

We revise th e upd ate algorithm for the smoothing 
len gth in Kawata _( 199£ ) who uses the algorithm suggested 
by Hernquist fc Katz_ ( 1989). The new al gorithm is based 

(|2000[). First, we count 



on that suggested by Thacker et al 



neighbour particles for the i-th particle with hi using the 
following smoothed kernel. 



W„n{r/h,) 



^W{4{r/h. 



if < r/h, < 3/4, 
3/4)) if 3/4 < r/h^ < 1, 



where r is the distance of the neighbour particle from the 
i-th particle, and W{x) is the same spline kernel as equa- 
tion (^) ^. The number of neighbours, A'^^ib.i, is no longer 
integer, but real number. This revision helps alleviate the 
discontinuity in the number of neighbours. 

Next, to avoid a sudden c hange in t h e smo othing length, 
we modify equation (14) of Kawata (1999) which deter- 



mines the evolution of the sm oothing l e ngth. Setting s = 



as. 



,1/3 



equation (14) of 



Kawata 



( 199£ ) is expressed 



h1{l — a -I- as) 



(9) 



where Ns — 40 is the desired number of neighbours, hi 



means the smoot hing len 



Kawata 



jthof the i-th particle at n step. 
(1999) corresponds to the case of 



Equation (14) of 
a = 0.5. Here we change a as a function of s 

0.2(1 + 3^) if s < 1, 
0.2(1 + l/s3) ifs>l. 



(10) 



In addition, 



Thacker et al. 



(2000) 



a modified 

Courant condition depending on the relative velocities to 
neighbour particles. Because this condition requires a sig- 
nificant amount of additional computational cost, once star 
formation and kinetic feedback are involved, we a b ando n 
this condition. As shown in Fig. 1 of Thacker et al. ( 2000 ), 
the above two algorithms lead to a dramatical improvement 
in keeping the number of neighbour particles constant, and 
the modified Courant condition adds a marginal effect. 

We follow the evolution of the smoothing length of star 
particles until the end of simulation, because we relax the 
instantaneous recycling for SNe II and consider the feedback 
and chemical enrichment from SNe la and intermediate mass 
stars. Even in the above new algorithm, involving kinetic 
feedback enables gas particles to escape from the position of 
star particles more quickly than the speed of the increase in 
the smoothing length of the star particles. In serious cases, 
it happens that there is no neighbour particle, and the mass 

^ Note that our s r nooth ing definition is different from that used 



in Thacker et al. 



twice as large as theirs. 



( |200o| ), in other words our smoothing length is 
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Figure 3. Comparison of the CMRs for the si mulation end- 
products and C oma Cluster galaxies (small dots, Bower, Lucey 
Ellis 1992a) in the aperture of 5 kpc. The circles connected 



by solid ines indicate the CMRs for the best model. The square 
(sta r ) de notes the result of model noIa-LM (illbla-HM) in Section 
O (O). The d ashed line shows the CMR fit ted to the Coma 



Cluster galaxies (Bower, Lucey, & Ellis L992b) 



conservation is broken due to the missing of the mass depo- 
sition from stars to gas. To solve this problem, we introduce 
the iteration process using equation once the number of 
neighbour particles is less than ten. Although this iteration 
requires an extra computational cost, it ensures complete 
mass conservation. 



3 APPLICATION TO ELLIPTICAL GALAXY 
FORMATION 

To see how the new code works, in this section, we ap- 
ply GCD+ to simulations of elliptical galaxy formation in 
the semi-cosmological model described in KOlb. KOlb ex- 
amined the mass-dependent properties of elliptical galaxies, 
such as the CMR, the Kormendy relation, and the [Mg/Fe]- 
magnitude relation. Our new code relaxes some of the sim- 
ple assumptions in KOlb, and becomes a more sophisticated 
code. Hence, the re-examination of those properties is an in- 
teresting application as well as an important test for GCD+. 



Figure 4. The CMRs for the simulation end-products in the 
99 kpc aperture. The circles connected by solid lines indicate the 
results of the best model. The square (star) denotes the result of 
model noIa-LM (illbla-HM). The dashed lines are the same as 
those in Fig. |3| 



3.1 Elliptical Galaxy Formation Model 

Following KOlb, we consider an isolated sphere, as a seed 
galaxy, on which small-scale density fluctuations corre- 
sponding to a cold dark matter (CDM) power spectrum are 
superi mposed. Here, we u se Bertschinger's software COS- 



MICS ( Bertschinger 1995) to generate initial density fluctu- 



ations. To incorporate the effects of fluctuations with longer 
wavelengths, the density of the sphere has been enhanced 
and a rigid rotation corresponding to a spin parameter. A, 
has been added. The initial conditions of this model are de- 
termined by the following four parameters: A, Aftot, erg, in, 
and Zc- The spin parameter is defined by 



A 



J\E 



1/2 



5/2 



(11) 



where J is the total angular momentum of the system, E is 
the total energy, and A/tot is the total mass of the sphere, 
which is composed of dark matter and gas; erg, in is the rms 
mass fluctuation in a sphere of radius 8 Mpc, which nor- 
malizes the amplitude of the CDM power spectrum; is the 
expected collapse redshift. If the top-hat density perturba- 
tion has an amplitude of 5i at the initial redshift, Zi, we ob- 
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tarn Zc 



1 9931) 
T99q) 



0.365,(1+2, 
rhus, when Zc is given 



1 approximately (e.g. | Padmanabhan 
5i at Zi is determined. Kawata 



is determined 
found that the seed galaxy which has a slow rotation 
corresponding to A = 0.02 and small-scale density fluctua- 
tions evolves into an elliptical-like system. Thus, we employ 
A = 0.02. This spin parameter is close to the minimum value 
possible in a CDM universe, according to the results of N 



body sim ulations (Barnes & Efstathiou 



1987 



Warren et al 



1992). We fix us, in = 0.5 and Zc — 3.5. Our simulations 
assume a flat universe (S7 = 1) with a baryon fraction of 
fib = 0.1 and a Hubble constant of Ho = 50 km s~^ Mpc~^. 
We carry out each simulation using 9171 particles for gas 
and dark matter, respectively. We simulate the evolution of 
each model from Zi = 40 to z = 0. 

The morphological evolution of all the models which 
are sho wn in this paper are similar to the evolution seen in 
Fig. 1 of Kawata (1999). A nearly spherical stellar system is 



formed at 2; = in all the models. Following KOlb, we focus 
on the chemical and photometric properties for the simu- 
lation end-products at z = 0. In our simulations the stel- 
lar particles each carry their own age and metallicity "tag" , 
which, when combined with population synthesis, enables us 
to derive the photometric properties of the simulated stel- 
lar systems. The photometric properties are derived by the 
same analysis as KOlb. De t ails of this analysis are described 
in Section 3.2 of [Kawata ( 2001a ) and Section 4 of KOlb. 

In this analysis, the spectral energy distribution (SED) 
of each stellar particle is assumed to be that of a single stellar 
population (SSP) that means a coeval and chemically homo- 
geneous assembly of stars. Since the observational data with 
which our results should be compared provide the luminos- 
ity distribution projected to a plane, we have to derive the 
projected distribution of SED from the three dimensional 
distribution of stellar particles. Finally, we obtain the x- y 
projected images as shown in Figure 5 of Kav 
when z-axis is set to be the initial rotation axis. 



( pOOlaD , 
We con- 
firmed that the results do not depend on the direction of 
the projection, because the end-products are nearly spher- 
ical. The flux of each stellar particle is smoothed using a 
Gaussian fllter with the fllter scale of 1/4 of the softening 
length of the stellar particle. These images provide similar 
information to the imaging data obtained in actual observa- 
tions. Thus, we can obtain various photometric properties 
from these images in the same way as in the analysis of ob- 
servational imaging data. We use the imag es similar to the 
one displayed in Figure 5 of Kawata (2001a), but employing 
a 1001 X 1001 pixel mesh to span the squared region with 
100 kpc on a side. 

We adopt the data of SSPs of Kodama fc A rimoto 97 

Kodama 



model (Kodama 



1997 



Kodama & Arimoto 



1997) 



& Arimoto 97 model supplies the database of SSPs with 
two types of IMF: {x, Mu, Mi) = (1.35, 60, 0.1) and (x 
Mu, Ml) = (1.1, 60, 0.1) in the definition of Section 
We adopt the data of SSPs with the IMF of (x, Mu, Mi) 
= (1.35, 60, 0.1), while we use the IMF of {x, Mu, M) = 
(1.35, 60, 0.2) in the numerical simulations. As mentioned 
in Section 4 of KOlb, this inconsistency does not affect the 
final photometric properties. 

The global photometric properties, such as the total 
luminosity, the colours, and the effective radius, are obtained 
from the projected image data. Following KOlb, the total 
magnitude and the effective radius are derived by fitting the 



O 

o o 



N 
N 



o 



o 

o 

-+ 
o 

CN 




9) ® — 



CD 

< 



o 



O) - 



□ 



-24 



-22 



-20 



Figure 5. The metallicitics (upper panel) and ages (lower panel) 
against the absolute V band magnitude for each model. The cir- 
cles connected by solid lines indicate the results of the best model. 
The square (star) denotes the result of model noIa-HM (illbla- 
LM). The open (filled) symbols denote the values evaluated in the 
5 kpc (99 kpc, spread over almost the whole galaxy) aperture. 



surface brightness profile to the r^^" law (eq. [11] of KOlb). 
We set the center to the position of a pixel which has the 
maximum V band luminosity. 

These procedures enable us to compare the properties 
for simulation end-products with those for observed ellipti- 
cal galaxies directly and quantitatively. In the following sec- 
tions, we begin by showing the best model which can repro- 
duce the CMR of observed elliptical galaxies, and examine 
what is responsible for the CMR based on the detailed anal- 
ysis of the simulation results. Next, we discuss the physical 
sizes of the simulation end-products, comparing them with 
the observed Kormendy relation. Finally, we examine the 
correlation between the total magnitude and the abundance 
ratio of Mg to Fe, i.e. the [Mg/Fe]-magnitude relation, which 
is sensitive to the star formation history. 



3.2 The Best Model: the Colour-Magnitude 
Relation 

As explained in Section ^, due to the lack of knowledge of 
the physics of star formation and SNe feedback, our code 
has a number of free parameters, such as the mass range 



8 D. Kawata and B.K. Gibson 




t (Gyr) 



t (Gyr) 



Figure 6. Time variation of tlie event rate of SNe II (thin l ines) and SNe la (thick lines) f or all the models. Points with error bars (at 
t = 13 Gyr) are taken from the observational SNe la rate by |C!apcllaro, Evans, fc Turatto 

I ( |l999| ). 



Table 1. Model Parameters 







Particle Mass (Mq) 


Softening (kpc) 






Name 


Mtot (Mq) 


Dark Matter 


Gas 


Dark Matter 


Gas 


SNe II 


SNe la 


HM 


4 X 10^2 


3.93 X 10** 


4.36 X 10'' 


4.28 


2.06 


yes 


yes 


MM 


8 X 10" 


7.85 X lO'^ 


8.72 X IQf' 


2.50 


1.20 


yes 


yes 


LM 


2 X 10" 


1.96 X 10^ 


2.18 X 10'' 


1.58 


0.758 


yes 


yes 


noIa-LM 


2 X 10" 


1.96 X 10^ 


2.18 X 10*5 


1.58 


0.758 


yes 


no 


illbla-LM 


2 X 10" 


1.96 X lO'^ 


2.18 X 10^ 


1.58 


0.758 


instantaneous 


1.5 Gyr delay 



of the IMF {Ml and Mu in Section 
parameters of /„ and esjv in Section 




and the feedback 
3.4 We fix these 
parameters so that the observed CMR is reproduced. To 
compare with the CMR, we simulate the evolution of seed 
galaxies with three different masses of Mtot = 4 x 10^^, 
8 X 10", and 2 x lO" Mq (hereafter, models HM, MM, and 
LM, respectively). The mass and spatial resolutions of each 
model are summarized in Table |l| 

The mass range of the IMF controls yields per mass 
of new born star particles. Higher or M\ leads to larger 
yields, and makes the metallicity of the end-products higher. 
We chose (Mu,Mi) = (60,0.2) Mq, which produces enough 
yields for the highest mass model (HM) to reproduce the 
observed red colour of galaxies with similar luminosities. 
As shown in KOlb, we also find that the stronger feedback 
causes the steeper slope in the CMR. The magnitude of the 
feedback effect is controlled by the parameters esjv and /„. 
We fix the parameter set of {esN,fv) = (0.1,0.002) such 
that the observed CMR is reproduced, when the above IMF 
is adopted. Accor ding to high resolution ID simulations of 
an SN remnant in Thornton et al. (199S), 90 % of the iin- 



tial SN energy is lost in radiation during its early expansion 
phase; this is not resolved in our simulations. The available 
SN ener gy of 10^° ergs which w e adopt is consistent with the 
result of Thornton et al. (199^ ), when a canonical initial SN 



energy of 10^ ergs is assumed. Thornton et al. (1998) also 
suggested that in the last stage of their simulation about 
~ 8.5 X lO''^ erg, i.e. 85 % of the available energy, is found 
in the form of kinetic energy. This corresponds to — 0.85 
in our definition, and is much larger than that we assumed 
here. However, the resolution limit of our simulation (a few 
kpc) is much larger than the physical scale of their simu- 
lation (~0.15 kpc). Unfortunately, little is know about the 
fraction of the kinetic energy which is still available to affect 
the interstell ar medium on su c h lar ge scales. Therefore, fv 
suggested by Thornton et al. (1995) is considered to be an 



upper limit, which does not invalidate the lower value which 
we adopted here. Although this parameter set is not a unique 
solution, it is one which best explains the observed CMR. 
Therefore, the set of models HM, MM, and LM with the 
above parameters is called the "best model" in this paper. 

Fig. ^ shows the CMR of the best model which repro- 
duces the observed CMR for galaxies in the Coma Cluster 
very well. The data for galaxies in the Coma Clus ter are 
the observed CMR of [Bower, Lucey, fc Ellis | (|l992^. Since 



there is no difference between SOs and ellipticals in the scal- 
ing relations which we now discuss, we do not dis t inguish 



between SOs and ellipticals. Bower, Lucey, fc Ellis ( 1992c ) 
supply U ~ V and V — K colours which refer to an aper- 
ture size of 11 arcsec, and the V band total magnitude de- 
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Table 2. Global Photometric Properties: Surface Brightness Profiles and Gradients 



V band K band 





Mv 






Mk 






[MB - m 


[Alog(Z/Z0)] 


Name 


(mag) 


n 


(kpc) 


(mag) 


n 


(kpc) 


1 A log r 


/ A log r 


HM 


-21.96 


5.80 


15.0 


-24.94 


6.51 


11.9 


-0.17 


-0.38 


MM 


-20.00 


3.63 


5.69 


-22.94 


3.86 


4.84 


-0.14 


-0.30 


LM 


-18.04 


1.80 


3.33 


-20.88 


1.71 


3.08 


0.00 


-0.01 


noIa-LM 


-18.58 


3.01 


1.32 


-21.62 


2.86 


1.09 


0.10 


-0.79 


illbla-LM 


-18.88 


7.26 


2.34 


-21.90 


7.47 


1.64 


-0.12 


-0.32 



Table 3. Photometric Properties and Stellar populations within Apertures 







D < 5 kpc 






D < 9 


3 kpc 




D <r^ 


s/2 










Age 








Age 






Name 


U-V 


V -K 


[Z/H] 


(Gyr) 


U-V 


V -K 


[Z/H] 


(Gyr) 


[Mg/Fc] 


[E/Fe 


HM 


1.56 


3.24 


0.17 


12.1 


1.29 


3.01 


-0.08 


12.1 


0.06 


0.10 


MM 


1.42 


3.13 


0.01 


12.1 


1.27 


2.97 


-0.14 


12.1 


0.10 


0.13 


LM 


1.25 


2.91 


-0.29 


12.2 


1.18 


2.84 


-0.34 


12.1 


0.14 


0.15 


noIa-LM 


1.32 


3.11 


0.25 


8.99 


1.25 


3.02 


0.13 


9.94 


0.05 


0.04 


illbla-LM 


1.47 


3.16 


0.08 


12.0 


1.36 


3.07 


-0.03 


12.0 


0.08 


0.11 



to 
d 



o 



d 




a small size. Observed elliptical galaxies (especially in lumi- 
nous galaxies) have colour gradients which make the colour 



12 



Log(M 



Figure 7. The ratio of the ejected gas mass (Mg Qw) at 2: = to 
the initial gas mass (Mg ;„;). The circles connected by solid lines 
indicate the results of the best model. The square (star) denotes 
the result of model noIa-LM (illbla-LM). Here, the ejected gas 
is defined as the sum of all the gas particles whose galactocentric 
radius is greater than 20 kpc. 



rived from a combination of their data and the literature. 
Throughout this paper we adopt the distance modulus of 
the Coma cluster of m — M = 34.7 mag; the Virgo distance 
modulus of m - M = 31.01 ( [Graham et al. |[l999| ) and the 
relative distance modul us of the Coma with respect to the 
Virgo of m - M = 3.69 ( [Bower, Lucey, fc Ellis | [l992b[ ) . This 
gives the luminosity distance of 87.1 Mpc for the Coma. We 
assume that the angular diameter distance equals the lu- 
minosity distance, because the redshift of the Coma cluster 
(z ^ 0.023) is nearly zero cosmologically. Then the aperture 
size of 11 arcsec at the distance of the Coma Cluster corre- 
sponds to ~ 5 kpc. The colours for simulation end-products 
are derived in the same aperture size as that used in the 
observations. 

KOlb claimed that the aperture effect should not be ig- 
nored, when discussing the CMR observed in an aperture of 



at the centre 


redder 


Peletier et al. 


199C|) 



The end-products in our simulations 
for modes HM and MM also have significant colour gradi- 
ents as seen in Table ^ Since colours within a fixed aperture 
yield colours in a more central region for larger galaxies, it 
is possible that the colours of large galaxies become redder 
than those of small galaxies, even if the mean colour of the 
whole galaxy is the same between the large and small galax- 
ies. Fig. ^ shows the CMR for the mean colour within the 
aperture of 99 kpc, which covers almost the entire galaxy in 
all the models which we have examined. As expected, the 
larger galaxy's colours are changed more dramatically by 
the aperture, which makes the slope between models HM 
and MM shallower than that between models MM and LM. 
The aperture effect contributes significantly to the slope of 
the CMR in Fig. [j] in the mass range b etween models HM 



1998) 



and MM (see also Kodama et al. 

As seen in KOlb, Table ^ shows that the colour and 
metallicity gradients become shallower, and Sersic law in- 
dex n (see eq. [11] of KOlb) in the surface brightness profile 
becomes smaller, with decreasing mass of the system. It may 
appear strange that the mean colour within the 99 kpc aper- 
ture is bluer than that within the 5 kpc aperture in model 
LM, although the colour gradient for model LM is flat. This 
is because Table ^ shows the gradients within the radius 
where the B-band magnitude is brighter than fiB = 24.5 
mag arcsec"^, following KOlb, and the gradient at larger 
radii becomes negative. 

The sequences of colours and line strengths among ellip- 
tical galaxies can almost equally well be attributed to ei ther 
age difference or metallicity differences (Worthey 1994). In 



our simulation this degeneracy can be broken completely, 
because we can directly analyse the metallicity and age in 
the simulation end-products. Fig. ^ shows the metallicity 
and age for each model. Here, the metallicities and ages are 
luminosity weighted values. Although the ages are all very 
old irrespective of the models and the apertures, the metal- 
licities exhibit similar behavior to the colours as a function 
of luminosity. Thus, we conclude that the slope of the CMR 
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Figure 8. The density-temperature distribution of the gas particles with higher (left), intermediate (middle), and lower (right) metallicity 
a,t z = 3.27 for models HM (upper panels) and LM (lower panels). The solid curves separate the region where the cooling times are 
shorter (upper region) and longer than the age of the universe a,t z = 3.27. 



of the best model in Fig. ^ is caused solely by the changes 
in metallicity. 

Fig. 1^ shows the histories of SNe II and SNe la event 
rates normalized by the total mass of the system. In our sim- 
ulations, the histories of SNe II events roughly trace those 
of star formation, because the lifetime of SNe II progenitors 
is shorter than the size of the bins of the histograms. We 
can clearly see that star formation ceases abruptly around 
t = 1.7 Gyr {z = 1.7), regardless of the mass of the system. 
We confirmed that gas particles begin to blow out from the 
stellar system, i.e. the galactic wind occurs when star for- 
mation stops. Subsequently, all the gas particles overcome 
the binding energy of the dark matter and stars, and escape 
from the system, as shown in Fig. 9 of KOlb. 

Fig. 1^ shows the mass fraction of the ejected gas as a 
function of the mass of the system. The ejected gas mass in 
Fig. ^ is defined as the total mass of all the gas particles 
whose galactocentric radius is greater than 20 kpc at z = 0. 
Comparing between models MM and LM, the lower mass 
system ejects a greater fraction of gas. Since the ejected gas 
cannot contribute to further metal enrichment, metal enrich- 
ment is more strongly suppressed in the lower mass system, 
leading to its lower metallicity in Fig. ^. On the other hand, 
between models HM and MM, the lower mass system ejects a 
slightly smaller fraction of the gas. In this luminosity range, 
the aperture effect appears to be a dominant factor in ex- 
plaining the slope of the mass-metallicity relation. However, 
there is a significant slope between models HM and MM in 
mean metallicity within the large aperture (Fig. ^) . We find 
that this is because there is another factor which contributes 
to the mass-metallicity relation. Fig. ^ plots density against 
temperature for the gas particles aX z = 3.27 for models HM 
and LM. At z = 3.27 the system is almost relaxed, the gas 



particles which stay in the lower right region in the panel 
are hot gas which has not yet cooled since they virialised. 
On the other hand, the gas particles which stay in the upper 
left (logpgas > —28 and logT* < 4) region represent the cold 
gas. The gas particles in the lower left region are escaping 
from the system. Comparing between models, the high mass 
system has more hot gas and less cold gas than the low mass 
system. The solid lines in the panels indicate where the cool- 
ing time equals the age of the universe at z — 3.27. Cooling is 
efficient in the upper region from this line, and inefficient in 
the lower region from this line. In the high mass system, the 
hot gas particles fall into where cooling is inefficient. There- 
fore, a large fraction of the gas cannot cool or infall into the 
central region where stars are forming, thus the hot gas is 
not affected by metal enrichment from stars (Fig. Since 
the amount of gas is small in the central region, relative to 
the amount of metals returned from stars, the metallicity 
of the gas in the central region increases rapidly. In other 
words, a low infall rate of the gas induced by inefficient cool- 
ing leads to efficient chemical enrichment. This is a similar 
mechanism to the infall model which was proposed to solve 
the G-dwarf problem of the clo sed box chem ical evolution 
model for the solar vicinity (e.g. Pagel 1997). On the other 
hand, in the low mass system the hot gas stays close to the 
line in Fig. ^. A large fraction of this gas which cools effi- 
ciently infalls into the central region, and suffers from metal 
enrichment. Since there is a large amount of gas in the star 
forming region, the metallicity of the gas in the central re- 
gion slowly increases. In other words, a high infall rate of gas, 
caused by efficient cooling, leads to inefficient chemical en- 
richment. Consequently, the high mass system makes metal 
rich stars more quickly than the low mass system (Fig. 
and the mean metallicity of the high mass system becomes 
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Figure 9. Metallicities against age for the star particles in the 
final stellar systems of models HM (upper panel) and LM (lower 
panel) . 



higher than that of the low mass system. Hence, we con- 
clude that not only the galactic wind, which ejects a larger 
fraction of gas in the lower mass system, but also the mass 
dependence of the gas infall rate, contributes to the slope in 
the mass-metallicity relation (Fig. |^ and the CMR (Fig. |^. 



3.3 The Role of SNe la 

To assess the relative importance of SNe la to SNe II during 
the evolution of elliptical galaxies, we have run an additional 
model which is similar to model LM, except that SNe la 
are not included. Hereafter, this model is called "noIa-LM". 
Figs. also show results of model noIa-LM. Comparison 
between models LM and noIa-LM demonstrates that SNe 
la are crucial sources for the galactic wind. For example, 
Fig. ^ shows that model noIa-LM does not have a clear 
cessation of star formation. Consequently, model noIa-LM 
ejects a smaller fraction of gas from the system (Fig. ^ , and 
has a higher metallicity and a redder colour than model LM 
(Figs. ^ and ^). In the U — V and Mv diagram, model nola- 
LM appears to follow the CMR of the best model. This is 
due to the younger age of model noIa-LM as seen in the lower 
panel of Fig. Model noIa-LM adopts the same feedback 
parameters as those of model LM. In addition, model nola- 
LM has a larger number of SNe than model LM, because 



more stars have formed. Thus, this result indicates that the 
galactic wind is caused mainly by SNe la. It is also notable 
that the event rate of SNe la exceeds that of SNe II, when 
the galactic wind occurs in models HM, MM, and LM. We 
therefore conclude that SNe la play a crucial role in driving 
and maintaining the galactic wind, and thus in contribut- 
ing to the evolution of elliptical galaxies. This conclusion 
is consistent with KOlb's suggestion, and our sophisticated 
treatment of SNe II and SNe la reinforces this idea. 



3.4 Instantaneous recycling SNe II and burst SNe 
la model 

KOlb adopted instantaneous recycling for SNe II and a sim- 
ple burst SNe la model, in which all the SNe la occur si- 
multaneously 1.5 Gyr after star formation. To compare the 
feedback effect of this simple description of SNe in KOlb with 
that of the more realistic SNe model used in this paper, we 
now examine another model which employs instantaneous 
recycling SNe II and 1.5 Gyr delay burst SNe la. The to- 
tal number of SNe II is given by the integral of number of 
stars with the mass higher than 8 Mq. The total number of 
SNe la expected in each star particle is assumed to be given 
by the sum of the MS-I-WD and RG-I-WD binary systems, 
which corresponds to the integral of equation (^) from t = 
to t = oo. We also assume that SNe la occur only in star par- 
ticles with log Z/Zq > -1.1, following the KTNOO model. 
Using this SNe model, we carry out the simulation with the 
same total mass as model LM. Hereafter, this model is called 
"illbla-LM". 

Fig. 1^ shows that the galactic wind ceases star forma- 
tion in this model, although the amount of ejected gas is 
much smaller than that for model LM (Fig. ^). As a re- 
sult, model illbla-LM has a higher metallicity and a redder 
colour than model LM (Figs. ^and|^. In both models LM 
and illbla-LM, the galactic wind occurs just after SNe la 
ignite. However, SNe la in model illbla-LM begin to occur 
later than those in model LM, due to the simple 1.5 Gyr 
delay, and model illbla-LM has a longer duration of star 
formation than model LM. As a result, in model illbla-LM 
the chemical enrichment progresses for a longer time, which 
leads to the higher metallicity and redder colour. Thus, the 
final properties of elliptical galaxies are sensitive to the time 
delay of SNe la, which is the main trigger of the galactic 
wind. It is also worth noting that the peak rate of star for- 
mation in model illbla-LM is significantly higher than that 
in model LM (Fig. This means that the suppression of 
star formation by continuous SNe II feedback is stronger 
than that from instantaneous recycling SNe II. These results 
demonstrate that we should not ignore the lifetime of both 
SNe II and SNe la progenitors, in studying the dynamical 
evolution of elliptical galaxies. 

3.5 The Kormendy Relation 

We now examine the physical size of the simulation end- 
products for all the models presented above. Fig. ^ shows 
a comparison of the Kormendy relations for the simulation 
end-products and the Coma Cluster galaxies both in the V 
and K ban ds. We refer to the data of the Coma Cluster 
galaxies of Pahre (199!:). In deriving the absolute magni- 



tude and effective radius in the kpc from the data set in 
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son of the theoretical model with the observed absorption 
line features. We compare our results with the observational 
results in terms of [Mg/Fe] rather than the absorption line 
features, because of convenience, although we plan to anal- 
yse the absorption line features from simulation results and 
compare them with the observational data more directly in 
future work. The observational results indicate that [Mg/Fe] 
correlates with the luminosity, i.e. the [Mg/Fe]-magnitude 
relation. Because Mg and Fe are mostly produced by SNe II 
and SNe la, respectively, and because SNe la have a longer 
delay than SNe II after the formation of stars, this corre- 
lation provides a strong co nstraint on the star formation 
histor v of elliptical galaxies ( Matteucci^ Ponzone, fc Gibson 
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Figure 10. Comparison of the Kormendy relations for the simu- 
lation en d-products and Coma Cluster galaxies (small dots, Pahrc 
"l99 S|) i i the V (upper panel) and K (lower panel) bands. The 



circles connected by solid lines indicate the Kormendy relation 
for the best model. The square (star) denotes the result of model 
noIa-LM (illbla-HM). 



Pahre 



(199£), we assume the same distance modulus as 



mentioned in Section 3.2, The best model shows a tendency 
for higher mass galaxies to have larger effective radii, which 
is qualitatively consistent with the observed Kormendy re- 
lation. However, the effective radii are systematically larger 
and the slope is slightly shallower than the observational 
data. In comparison among the low mass systems, model 
LM shows a significantly larger effective radius than models 
noIa-LM and illbla-LM. Sections 3. J and 3.4 show that the 



effect of SNe on star formation in the models decreases as: 
LM, illbla-LM, and noIa-LM, and the effective radius de- 
creases in the same order. In other words, stronger feedback 
causes an expansion of the system. Due to this expansion ef- 
fect, the best model fails to explain the observed Kormendy 
relation, which is the same probl em as that found in KOlb 
(see also Chiosi fc Carraro 2002). 



3.6 The [Mg/Fe]-Magnitude Relation 

The abundance ratio of [Mg/Fe] of el liptical galaxies has 



been well-studied by several groups fVVorthey, Faber, & 


Gonzalez 1992 ^ ]0rgcnsen 1999; 


Kuntschner ^000|; Irager 


et al. 2000lj|; [Tcrlcvich & Forbes 


|2002D. Strictly speaking. 



[Mg/Fc] is not "observed", but is derived by the compari- 



199S 



2002 



Thomas, Grcggio, fc Bender 



Tantalo fc Chiosi 2002). In Fig, 



: compare [E/Fe] 

for the sim ulation end-pro d ucts o f all the models with those 



I199C 

yTw 



Romano et al. 



, . .|2000a|) 

ba nd magnitude s in tTrager et al.] (2000b). In both the data 



derived by Trager et al. (2000a|). We use the absolute B 



of Trager et al. (2000a) and the simulation results, [E/Fe 
i s meas ured in the re,B /'i aperture 

(^00040) ad opt [E/Fe] instead of [Mg/Fe] (see also [Tantalo 



et al. 



998ej ), we derived [E/Fe] as well as [Mg/Fe] |]. The 



definition of [E/Fe] is as follows: 
[E/Fe] = ([Z/H] - [Fe/H])/A, 



(12) 



where A = 0.929 for model 4 of Trager et al 



we refer. In model 4 of Trager et al. (|2000a), [E/Fe] indicates 



(2000a) which 



the mean enhancement of "E" group, which contains C, O, 
Ne, Na, Mg, Si, and S, with respective to the Fe-peak ele- 
ments. The quantity of [E/Fe] should be similar to [Mg/Fe], 
because most of the E group elements are nucleosyntheti- 
cally linked. Table |^ summarises luminosity weighted [E/Fe] 
and [Mg/Fe] for all the model, and shows that the differ- 
ence between [E/Fe] and [Mg/Fe] is less than 0.04. There- 
fore, we consider that it is safe to use the [E/Fe]-magnitude 
relation i n order to ex a mine t he [Mg/Fe]-magnitude rela- 
tion. The Trager et al. (2000a) data show a clear tendency 
that [E/Fc] increases with the galactic luminosity. However, 
the best model is incapable of reproducing this tendency 
and has almost constant [E/Fe] regardless of the luminosity. 
Models noIa-LM and illbla-LM have similar [E/Fe] to the 
best model. As mentioned above, these [E/Fe] abundance 
ratios reflect the star formation history shown in Fig. ^. In 
all the models except noIaLM, the galactic wind causes the 
cessation of star formation, and star formation stops before 
the metal enrichment by SNe la progresses. Since [E/Fe] is 
determined only by SNe II yields, these models have simi- 
lar [E/Fe] to that for model noIa-LM. Finally, it is revealed 
that the best model to reproduce the observed CMR leads 
to a short period of star formation and constant [E/Fe], ir- 
respective of the mass of the system, which is inconsistent 
with the [E/Fe]-magnitude relation derived from observa- 
tions. Because the results of [E/Fe] are almost identical to 
those of [Mg/Fe], we focus on [Mg/Fe] in the following dis- 
cussion. 

A more detailed look shows that the best model has 
a slight slope which is opposite to the [Mg/Fe]-magnitude 
relation derived from observational data, and model illbla- 



Noto that Fig. 12 of KOlb compares [M 



end-products with [E/Fe] of Trager et al. 



; /Fe] of the simulation 



( 2000a| ). 
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LM (noIa-LM) has a little lower [Mg/Fe] than model LM 
(illbla-LM). These trends also reflect their star formation 
histories. Fig. ^ shows [Mg/Fe] vs [Fe/H] of all star particles 
for all the models. In all the models, there is large scatter in 
[Mg/Fe] at [Fe/H]< —1, and the mean [Mg/Fe] approaches 
0.1 at [Fe/H]> -1 . Model illbla-LM has the least scatter, 
compared to the other models, because of the assumption of 
instantaneous recycling for SNe II. To interpret Fig. the 
mean [Mg/Fe] for the total SNe II yields is useful, which is 
calculated by 



([Mg/Fe]„(Z)) 



log 



Jg° MMg,ii(m, Z)^{m)dm 
° MFc,ii(m, Z)^{m)dm 



-\og{Mg/Fe)Q. 



(13) 



This value equals [Mg/Fe] of the SNe II yields in 
model illbla-LM, and depends on the metallicity, such as 
([Mg/Fe]„(lO-*Z0)) ~ 0.3, ([Mg/Fe]ii(lO-2Z0)) ~ 0.08, 
and ([Mg/Fe]ii(Z0)> ~ 0.1. Therefore, the [Mg/Fe] ver- 
sus [Fe/H] distribution in model illbla-LM can be eas- 
ily explained as follows. Some low metallicity stars with 
0.3 > [Mg/Fe] > 0.1 formed from the gas enriched by the ex- 
tremely low metallicity {Z ~ 1O~''Z0) stars. The other stars 
( [Mg/Fe] ~ 0.1) are enriched mainly by stars with Z > 
n. The other models follow the same trend as model illbla- 
LM, although the y have much la rger scatter in [Mg/Fe]. As 
see n in Fig. 3 of jcibson | (|l997| ), [Mg/Fe] of SNe II yields 



by Woosley & Weaver 



( [1995| ) also depends on the progeni- 
tor mass (e.g. for solar metallicity stars, [Mg/Fe]=1.2 at 40 
M0 and [Mg/Fe]= -0.59 at 11 M0). Although the differ- 
ence in the life-times of SNe II progenitors is small, this time 
difference causes the inhomogeneous enrichment for the in- 
terstellar medium, and makes the scatter seen in Fig. [l2| 
Some fraction of stars with [Mg/Fe] < and [Fe/H]> —1 
are enriched by SNe la whose yields have [Mg/Fe]= —1.5. 
Model HM appears to have larger fraction of those stars than 
model LM. In addition, model LM stops star formation be- 
fore enrichment by SNe II whose yields have [Mg/Fe] ~ 0.1 
becomes dominant. These two effects lead to the slope of the 
best model seen in Fig. |ll|. 

The reason why [Mg/Fe] of model noIa-LM is smaller 
than that of model illbla-LM is a little more complicated 
- it is due to the enrichment from intermediate mass stars. 
In theory, with the exception of carbon, nitrogen, and oxy- 
gen, the abundance pattern of the yields for intermediate 
mass stars is the same as the abundance pattern when stars 
are born, i.e. the initial abundance pattern. However, the 
abundance pattern for intermediate mass stars is set to be 
the solar abundance patte rn in our code, due to the reason 
mentioned in Section ^.3.2 . As time progresses, the yields for 
intermediate mass stars becomes more important, because 
more stars with high metallicities finish the main-sequence 
phase. Therefore, model noIa-LM which has the longest 
duration of star formation makes a significant amount of 
stars which suffer from the yields with the "artificial" so- 
lar abundance pattern. These stars are seen in Fig. ^ as 
stars with high metallicity ([Fe/H]> 0.2) and slightly lower 
[Mg/Fe] (around 0.05). As a result, model noIa-LM has lower 



5 The stars with [Mg/Fe]= -0.18 and [Fe/H]< -1 arc enriched 
only by zero metallicity stars, because of ( [Mg/Fe] ii(0)) ~ —0.18. 
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Figure 11. Comparison of the simulation end- products with 
observed early-type galaxies (Trager et al. I I2OO0II ) in the [E/Fe] 
vs. Mb diagram. The circles connected by solid lines indicate the 
the best model. The square (star) denotes the result of model 
noIa-LM (illbla-HM). 



mean [Mg/Fc] than models LM and illbla-LM which have 
a shorter duration of star formation. Although inclusion of 
SNe la makes this effect negligible, this is caused by the nec- 
essarily simple recipe for the yields for intermediate mass 
stars, which should be improved in the near future. 



3.7 Discussion 

We have applied our new chemo-dynamical evolution code, 
GCD+, to simulations of elliptical galaxy formation, and con- 
firmed the results of KOlb for the CMR, the Kormendy re- 
lation, and the [Mg/Fe]-magnitude relation. Tables □ and 
^ summarise the photometric properties and stellar popu- 
lation for all the models presented. The most fundamental 
difference between our new code and the one used in KOlb is 
its taking into account of the lifetime of stars and the adop- 
tion of the new SNe la model proposed by KTNOO. We have 
obtained a more realistic history of SNe la rate than KOlb 
(Fig. |). The history of the SNe la rate follows the SNe la 
progenitor model of KTNOO. In the "best model" where star 
formation is stopped by galactic winds at high redshift, the 
histories of the SNe la rate have two peaks. These two peaks 
are caused by SNe la from the MS-I-WD and RG-(-WD sys- 
tems (Section 2.3.3). Their evolution is similar to that shown 



in Fig. 5 of KTNOO. After the second peak, this model pre- 
dicts the continuous explosion of SNe la from the RG-I-WD 
system until z — 0. Capellaro, Evans, & Turatto (199t) 



provides the observed SNe rates at nearby galaxies. This 
SNe II rate for elliptical galaxies gives an upper limit, and 
is consistent with the best model where there is no SNe II. 
Points with error bars in Fig. ^ show their observed SNe la 
rates which are transfered from units of SNu to units used in 
Fig. phased on the B band luminosity for each model. In the 
best model SNe la rates are in good agreement with the ob- 
served rates. This result which comes from a self-consistent 
chemo-dynamical evolution provides strong support for the 
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Figure 12. The [Mg/Fe] against [Fe/H] of the star particles in the final stellar system for all the models. 



KTNOO model which is based on pure chemical evolution 
studies. 

Our sophisticated SNe la model has clarified that the 
feedback effect of SNe la plays a crucial role in the evolution 
of elliptical galaxies rather than that of SNe II. The lifetimes 
of progenitor stars of SNe II are much shorter than those of 
SNe la. As shown in Fig. ^, SNe II occur immediately (be- 
tween 4 and 40 Myr for stars with Z=0.02) after star forma- 
tion events. On the other hand, there is a significant delay 
(> 0.7 Gyr for stars with Z=0.02) between star formation 
and SNe la. Therefore, SNe II occur in dense gaseous envi- 
ronments still actively forming stars, whereas SNe la act on 
tenuous gas left after star formation. This difference makes 
SNe la the major trigger of galactic winds. Comparison be- 
tween models noIa-LM and LM shows this clearly in Section 



3.3. Also model illbla-LM demonstrates that the simplified 
model for SNe II and SNe la underestimates significantly the 
SNe effect on the dynamical evolution of elliptical galaxies. 

The difference in SNe II and la model between KOlb 
and this paper is one of the reasons why our new model 
requires only 10^'' ergs and — 0.002 as an SN feedback 
energy to reproduce the observed CMR, whereas KOlb re- 
quires 4 X 10^^ ergs and /i, — 0.9. This difference is caused 
by the pnmhinatinn nf snme physirpil prnrpsses, such ss the 



We find that mass dependences of both the galactic 
wind and the gas infall rate are responsible for the slope 
of the CMR. The conventional galactic wind scenario con- 
siders only the former mechanism. Our numerical simulation 
shows that the higher mass system has a higher virial tem- 
perature and lower efficiency of radiative cooling. Therefore, 
the gas infall rate decreases with increasing mass of the sys- 
tem. This lower gas infall rate leads to the increase in the 
fraction of the high metaUicity stars (Fig ^ . As a result, the 
mean metaUicity of the higher mass system becomes higher 
than that of the lower mass system. This is a novel mecha- 
nism contributing to the CMR. 

Also in the [Mg/Fe]-magnitude relation, we reach the 
same conclusion as KOlb qualitatively. However, the value 
of [Mg/Fe] is systematically different from that in KOlb. 
This is due to the difference in the stellar vields of SNe I I. 
We use the yields calculated by [Woosley fc Weaver _( 1995 ) , 



wherea s KOlb adopts thos e by [Notomo et al. J 1997 ). Al - 
though Woosley fc Weaver ( 1995 ) and Notomo et al 



(1997 



are two representative computations for SNe II nucleosyn- 
th esis, it is known that the a mount of iron yield derived 



Woosley fc Weaver 



that in 



Woosley fc Weaver 



IMF, ^NTc TT anrl ClNo Tg focHhgr-lf anrl tVio rr^r^lino- fiinr- 



Wooslcy, fc Weaver 



(1995) is significantly larger than 
Notomo et al. | ( 1997), and therefore the yields of 
(|1995 ) lead to smaller [Mg/Fe] (|Timmes. ' 



19951; jcibson ||l997|; [Gibson, Lowenstein, 



tion, which are updated in the new code. Coniparirip- the 



Mushtzky 1 11997| ; [Thomas, Greggio, fc Bender ||199^ ). The 



results, our new code requires a more reasonable amount of 
the SN feedback energy than that in KOlb, to reproduce the 
observed CMR. However, we might still ignore some phys- 
ical processes which contribute to the CMR, and be able 
to further reduce the feedback energy. For example, the UV 
background radiation is a possible candidate as such physi- 
cal process, as suggested by KOlb. A recent semi-analytical 
study demonstrates that the UV background rad iation has a 



similar eff ect to the SNe feedback on the CMR ( Nagashima 
& GouifrilioOll). 



theoretical stellar yields still have a large uncertainty, and 
to discuss the zero-points of [Mg/Fe] for simulation end- 
products we need a more reliable yield model. Nevertheless, 
we can discuss the observed slope of the [Mg/Fe]-magnitude 
relation, and our model so far cannot explain this. 

To explain the observed slope, it is a plausible scenario 
that low mass systems have a longer duration of star for- 
mation a nd are enriched by SNe la more than high mass 



systems (^aber, Worthey, fc G onzalez 



1992 



Mattcucci 



1994 ; Matteucci, Ponzone, fc Gibson |199^ ). KOlb suggests 
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that the UV background radiation is a possible candidate 
to reahse the scenario, because it suppresses coohng and 
star form ation more strongly in lower mass systems (Ef- 
stathioi^ 1992), and is expected to prolong the duration 
of star formation. Since we have developed a sophisticated 
chemo-dynamical evolution code, it is important to intro- 
duce the UV background radiation into GCD+, and examine 
the effect of the UV background radiation on galaxy evolu- 
tion. In a forthcoming paper, we hope to throw a light on 
this issue. 



4 CONCLUSION 

We have developed a new chemo-dynamical evolution code 
called GCD+. The code includes both SNe II and SNe la mod- 
eling taking into account the lifetime of progenitor stars, 
and the chemical enrichment from intermediate mass stars. 
In this paper we have applied the code to simulations of 
elliptical galaxy formation described as the evolution of an 
isolated seed galaxy. We have shown that GCD+ is a useful 
and unique tool which enables us to compare simulation re- 
sults with the observational data directly and quantitatively, 
and gives us new insights into the effects of dynamics on the 
chemical evolution of galaxies. Encouraged by the success of 
this study, we plan to apply GCD+ to high-resolution cosmo- 
logical simulations in the near future. In the real universe, 
even isolated galaxies at z = are expected to experience 
interactions wit h the another galaxies in the past. In fact, 
(1998) suggests that interactions between galax- 



Gnedin 



ies are a dominant mechanism to transport metals from the 
galaxies to intergalactic medium rather than SNe feedback. 
Therefore, using simulations covering a cosmological scale, 
it will be valuable to examine what physical process, or com- 
bination of processes, induces the tight observed CMR. In 
addition, cosmological simulations using GCD+ can provide 
detailed metallicity distributions not only in the stellar com- 
ponent of individual galaxies, but also in the intergalactic 
medium, such as Lyman-a clouds and th e ho t gas in clus- 
ters of galaxies. As mentioned in Section |2.3| , the different 
chemical elements have progenitors with different masses, 
and different mass stars have different lifetimes. It is highly 
probable that the abundance ratios are fairly sensitive to 
the metal transport mechanism and its epoch. New obser- 
vational facilities will offer a wealth of information about the 
distribution of metallicity, as well as the abundance ratio of 
each element, such as [N/O], [Mg/Fe] and [Si/Fe], with un- 
precedented accuracy. Studies comparing new observations 
with our future cosmological simulation results will show 
clearly which metal transport mechanism is dominant, and 
thus contribute to understanding the formation and evolu- 
tion of galaxies as well as clusters of galaxies from a new 
perspective. 
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